#include "ani.h"
#define namelength      100   /* number of characters max. in species name */
#define down            2
#define over            60

#define ibmpc0          false
#define ansi0           true
#define vt520           false

typedef double *vector;
typedef long *intvector;
typedef Char naym[namelength];
typedef long longer[6];

typedef struct node {
  struct node *next, *back;
  boolean tip;
  long number;
  naym nayme;
  double v;
  long xcoord, ycoord, ymin, ymax;
} node;

typedef struct tree {
  node **nodep;
  node *start;
} tree;


Static FILE *infile, *outfile, *treefile;
Static long numsp, numsp1, numsp2, inseed, outgrno, col, datasets, ith, j;
Static vector *x;
Static intvector *reps;
Static naym *names;
Static boolean jumble, lower, upper, outgropt, replicates, trout,
               printdata, progress, treeprint, mulsets, ibmpc, vt52, ansi,
               njoin;
Static tree curtree;
Static longer seed;
Static long *enterorder;

/* Local variables for maketree, propagated globally for C version: */
node **cluster;


openfile(fp,filename,mode,application,perm)
FILE **fp;
char *filename;
char *mode;
char *application;
char *perm;
{
  FILE *of;
  char file[100];
  strcpy(file,filename);
  while (1){
    of = fopen(file,mode);
    if (of)
      break;
    else {
      switch (*mode){
      case 'r':
        printf("%s:  can't read %s\n",application,file);
	file[0] = '\0';
        while (file[0] =='\0'){
          printf("Please enter a new filename>");
          gets(file);}
        break;
      case 'w':
      case 'a':
        printf("%s: can't write %s\n",application,file);
	file[0] = '\0';
        while (file[0] =='\0'){
          printf("Please enter a new filename>");
          gets(file);}
        break;
      }
    }
  }
  *fp=of;
  if (perm != NULL)
    strcpy(perm,file);
}



double randum(seed)
long *seed;
{
  /* random number generator -- slow but machine independent */
  long i, j, k, sum;
  longer mult, newseed;
  double x;

  mult[0] = 13;
  mult[1] = 24;
  mult[2] = 22;
  mult[3] = 6;
  for (i = 0; i <= 5; i++)
    newseed[i] = 0;
  for (i = 0; i <= 5; i++) {
    sum = newseed[i];
    k = i;
    if (i > 3)
      k = 3;
    for (j = 0; j <= k; j++)
      sum += mult[j] * seed[i - j];
    newseed[i] = sum;
    for (j = i; j <= 4; j++) {
      newseed[j + 1] += newseed[j] / 64;
      newseed[j] &= 63;
    }
  }
  memcpy(seed, newseed, sizeof(longer));
  seed[5] &= 3;
  x = 0.0;
  for (i = 0; i <= 5; i++)
    x = x / 64.0 + seed[i];
  x /= 4.0;
  return x;
}  /* randum */




void uppercase(ch)
Char *ch;
{
    *ch = isupper(*ch) ? (*ch) : toupper(*ch);
}  /* uppercase */


void getnums()
{
  /* read species numbers or number of bootstrap data sets */
  fscanf(infile, "%ld", &numsp);
  fprintf(outfile, "\n%4ld Populations\n", numsp);
  numsp1 = numsp + 1;
  numsp2 = numsp * 2 - 2;
}  /* getnums */

void getoptions()
{
  /* interactively set options */
  long i, inseed0;
  Char ch;
  boolean done1;

 // fprintf(outfile, "\nNeighbor-Joining/UPGMA method version %s\n\n",VERSION);
  putchar('\n');
  jumble = false;
  lower = false;
  outgrno = 1;
  outgropt = false;
  replicates = false;
  trout = true;
  upper = false;
  printdata = false;
  progress = true;
  treeprint = true;
  njoin = true;
  for(;;) {/*
    printf(ansi ? "\033[2J\033[H" :
           vt52 ? "\033E\033H"    : "\n");
    printf("\nNeighbor-Joining/UPGMA method version 3.5\n\n");
    printf("Settings for this run:\n");
    printf("  N       Neighbor-joining or UPGMA tree?  %s\n",
           (njoin ? "Neighbor-joining" : "UPGMA"));
    if (njoin) {
      printf("  O                        Outgroup root?");
      if (outgropt)
        printf("  Yes, at species number%3ld\n", outgrno);
      else
        printf("  No, use as outgroup species%3ld\n", outgrno);
    }
    printf("  L         Lower-triangular data matrix?  %s\n",
           (lower ? "Yes" : "No"));
    printf("  R         Upper-triangular data matrix?  %s\n",
           (upper ? "Yes" : "No"));
    printf("  S                        Subreplicates?  %s\n",
           (replicates ? "Yes" : "No"));
    printf("  J     Randomize input order of species?");
    if (jumble)
      printf("  Yes (random number seed =%8ld)\n", inseed0);
    else
      printf("  No. Use input order\n");
    printf("  M           Analyze multiple data sets?");
    if (mulsets)
      printf("  Yes, %2ld sets\n", datasets);
    else
      printf("  No\n");
    printf("  0   Terminal type (IBM PC, VT52, ANSI)?  %s\n",
           (ibmpc ? "IBM PC" : ansi ? "ANSI" : vt52 ? "VT52" : "(none"));

    printf("  1    Print out the data at start of run  %s\n",
           (printdata ? "Yes" : "No"));
    printf("  2  Print indications of progress of run  %s\n",
           (progress ? "Yes" : "No"));
    printf("  3                        Print out tree  %s\n",
           (treeprint ? "Yes" : "No"));
    printf("  4       Write out trees onto tree file?  %s\n",
           (trout ? "Yes" : "No"));
    printf( "\nAre these settings correct?");
    //printf(" (type Y or the letter for one to change)\n");
    //scanf("%c%*[^\n]", &ch);
    //getchar();
    //if (ch == '\n')
    //  ch = ' ';
    //uppercase(&ch);*/
	ch = 'Y';
    if  (ch == 'Y')
      {njoin = !njoin;break;}
    if (strchr("NJOULRSM01234",ch) != NULL){
      switch (ch) {
	
      case 'J':
	jumble = !jumble;
	if (jumble) {
	  do {
	    printf("Random number seed (must be odd)?\n");
	    scanf("%ld%*[^\n]", &inseed);
	    getchar();
	  } while (!(inseed & 1));
	  inseed0 = inseed;
	  for (i = 0; i <= 5; i++)
	    seed[i] = 0;
	  i = 0;
	  do {
	    seed[i] = inseed & 63;
	    inseed /= 64;
	    i++;
	  } while (inseed != 0);
	}
	break;
	
      case 'L':
	lower = !lower;
	break;
	
      case 'O':
	outgropt = !outgropt;
	if (outgropt) {
	  done1 = true;
	  do {
	    printf("Type number of the outgroup:\n");
	    scanf("%ld%*[^\n]", &outgrno);
	    getchar();
	    done1 = (outgrno >= 1 && outgrno <= numsp);
	    if (!done1) {
	      printf("BAD OUTGROUP NUMBER: %4ld\n", outgrno);
	      printf("  Must be in range 1 -%2ld\n", numsp);
	    }
	  } while (done1 != true);
	}
        else
          outgrno = 1;
	break;
	
      case 'R':
	upper = !upper;
	break;
	
      case 'S':
	replicates = !replicates;
	break;
	
      case 'N':
	njoin = !njoin;
	break;
	
      case 'M':
	mulsets = !mulsets;
	if (mulsets) {
	  done1 = false;
	  do {
	    printf("How many data sets?\n");
	    scanf("%ld%*[^\n]", &datasets);
	    getchar();
	    done1 = (datasets >= 1);
	    if (!done1)
	      printf("BAD DATA SETS NUMBER:  it must be greater than 1\n");
	  } while (done1 != true);
	}
	break;
	
      case '0':
	if (ibmpc) {
	  ibmpc = false;
	  vt52 = true;
	} else {
	  if (vt52) {
	    vt52 = false;
	    ansi = true;
	  } else if (ansi)
	    ansi = false;
	  else
	    ibmpc = true;
	}
	break;
	
      case '1':
	printdata = !printdata;
	break;
	
      case '2':
	progress = !progress;
	break;
	
      case '3':
	treeprint = !treeprint;
	break;
	
      case '4':
	trout = !trout;
	break;
      }
    } else
      printf("Not a possible option!\n");
  }
}  /* getoptions */


void doinit()
{
  /* initializes variables */
  getnums();
  getoptions();
}  /* doinit */


void inputoptions()
{
  /* read options information */
  Char ch;
  long cursp;

  if (ith != 1) {
    if (eoln(infile)) {
      fscanf(infile, "%*[^\n]");
      getc(infile);
    }
    fscanf(infile, "%ld", &cursp);
    if (cursp != numsp) {
      printf("\nERROR: INCONSISTENT NUMBER OF SPECIES IN DATA SET %4ld\n",
             ith);
      exit(-1);
    }
  }
  while (!eoln(infile)) {
    ch = getc(infile);
    if (ch == '\n')
      ch = ' ';
    uppercase(&ch);
    if (ch != ' ') {
      printf("BAD OPTION CHARACTER: %c\n", ch);
      exit(-1);
    }
  }
  putc('\n', outfile);
  if (njoin)
    fprintf(outfile, " Neighbor-joining method\n");
  else
    fprintf(outfile, " UPGMA method\n");
  //fprintf(outfile, "\n Negative branch lengths allowed\n\n");
}  /* inputoptions */


void getinput()
{
  /* reads the input data */
  inputoptions();
}  /* getinput */




void setuptree(a)
tree *a;
{
  /* initialize a tree */
  long i, j;
  node *p, *q;

  a->nodep = (node **)Malloc((numsp2+1)*sizeof(node *));
  for (i = 0; i < numsp; i++) {
    a->nodep[i] = (node *)Malloc(sizeof(node));
    a->nodep[i]->tip = true;
    a->nodep[i]->number = i + 1;
    memcpy(a->nodep[i]->nayme, names[i], sizeof(naym));
    a->nodep[i]->v = 0.0;
  }
  for (i = numsp1; i <= numsp2; i++) {
    q = NULL;
    for (j = 1; j <= 3; j++) {
      p = (node *)Malloc(sizeof(node));
      p->tip = false;
      p->number = i;
      p->next = q;
      q = p;
    }
    p->next->next->next = p;
    a->nodep[i - 1] = p;
  }
  q = NULL;
  for (j = 1; j <= 2; j++) {
    p = (node *)Malloc(sizeof(node));
    p->tip = false;
    p->number = numsp2 + 1;
    p->next = q;
    q = p;
  }
  p->next->next = p;
  a->nodep[numsp2] = p;
  a->start = a->nodep[0];
}  /* setuptree */

void getdata()
{
  /* read in distance matrix */
  long i, j, k, columns;
  boolean skipit, skipother;

  columns = 6;
  if (printdata) {
    fprintf(outfile, "\nName                       Distances");
    if (replicates)
      fprintf(outfile, " (replicates)");
    fprintf(outfile, "\n----                       ---------");
    if (replicates)
      fprintf(outfile, "-------------");
    fprintf(outfile, "\n\n");
  }
  for (i = 0; i < numsp; i++) {
    x[i][i] = 0.0;
    fscanf(infile, "%*[^\n]");
    getc(infile);
    /*
    for (j = 0; j < namelength; j++) {
      if (eoln(infile))
	names[i][j] = ' ';
      else {
	names[i][j] = getc(infile);
	if (names[i][j] == '\n')
	  names[i][j] = ' ';
      }
    }
    */
    fscanf(infile,"%s",names[i]);
    for (j = 0; j < numsp; j++) {
      skipit = ((lower && j + 1 >= i + 1) || (upper && j + 1 <= i + 1));
      skipother = ((lower && i + 1 >= j + 1 )|| (upper && i + 1 <= j + 1));
      if (!skipit) {
	if (eoln(infile)) {
	  fscanf(infile, "%*[^\n]");
	  getc(infile);
	}
	fscanf(infile, "%lf", &x[i][j]);
	if (replicates) {
	  if (eoln(infile)) {
	    fscanf(infile, "%*[^\n]");
	    getc(infile);
	  }
	  fscanf(infile, "%ld", &reps[i][j]);
	} else
	  reps[i][j] = 1;
      }
      if (!skipit && skipother) {
	x[j][i] = x[i][j];
	reps[j][i] = reps[i][j];
      }
    }
  }
  fscanf(infile, "%*[^\n]");
  getc(infile);
  if (!printdata)
    return;
  for (i = 0; i < numsp; i++) {
    /*
    for (j = 0; j < namelength; j++)
      putc(names[i][j], outfile);
      */
    fprintf(outfile,"%s",names[i]);
    putc(' ', outfile);
    for (j = 1; j <= numsp; j++) {
      fprintf(outfile, "%10.5f", x[i][j - 1]);
      if (j % columns == 0 && j < numsp) {
        putc('\n', outfile);
        for (k = 1; k <= namelength + 1; k++)
          putc(' ', outfile);
      }
    }
    putc('\n', outfile);
  }
  putc('\n', outfile);
}  /* getdata */

Local Void hookup(p, q)
node *p, *q;
{
  /* hook together two nodes */
  p->back = q;
  q->back = p;
}  /* hookup */


void coordinates(p, lengthsum,tipy,tipmax)
node *p;
double lengthsum;
long *tipy;
double *tipmax;
{
  /* establishes coordinates of nodes */
  node *q, *first, *last;

  if (p->tip) {
    p->xcoord = (long)(over * lengthsum + 0.5);
    p->ycoord = (*tipy);
    p->ymin = (*tipy);
    p->ymax = (*tipy);
    (*tipy) += down;
    if (lengthsum > (*tipmax))
      (*tipmax) = lengthsum;
    return;
  }
  q = p->next;
  do {
    coordinates(q->back, lengthsum + q->v, tipy,tipmax);
    q = q->next;
  } while ((p == curtree.start || p != q)
        && (p != curtree.start || p->next != q));
  first = p->next->back;
  q = p;
  while (q->next != p)
    q = q->next;
  last = q->back;
  p->xcoord = (long)(over * lengthsum + 0.5);
  if (p == curtree.start) {
    if (njoin)
      p->ycoord = p->next->next->back->ycoord;
    else
      p->ycoord = (p->back->ycoord + p->next->back->ycoord) / 2;
  } else
    p->ycoord = (first->ycoord + last->ycoord) / 2;
  p->ymin = first->ymin;
  p->ymax = last->ymax;
}  /* coordinates */

void drawline(i, scale)
long i;
double scale;
{
  /* draws one row of the tree diagram by moving up tree */
  node *p, *q;
  long n, j;
  boolean extra, trif;
  node *r, *first, *last;
  boolean done;

  p = curtree.start;
  q = curtree.start;
  extra = false;
  trif = false;
  if (i == p->ycoord && p == curtree.start) {
    if (p->number - numsp >= 10)
      fprintf(outfile, "-%2ld", p->number - numsp);
    else
      fprintf(outfile, "--%ld", p->number - numsp);
    extra = true;
    trif = true;
  } else
    fprintf(outfile, "  ");
  do {
    if (!p->tip) {
      r = p->next;
      done = false;
      do {
        if (i >= r->back->ymin && i <= r->back->ymax) {
          q = r->back;
          done = true;
        }
        r = r->next;
      } while (!(done || (p != curtree.start && r == p) ||
                 (p == curtree.start && r == p->next)));
      first = p->next->back;
      r = p;
      while (r->next != p)
        r = r->next;
      last = r->back;
      if (p == curtree.start)
        last = p->back;
    }
    done = (p->tip || p == q);
    n = (long)(scale * (q->xcoord - p->xcoord) + 0.5);
    if (n < 3 && !q->tip)
      n = 3;
    if (extra) {
      n--;
      extra = false;
    }
    if (q->ycoord == i && !done) {
      if (p->ycoord != q->ycoord)
        putc('+', outfile);
      if (trif) {
        n++;
        trif = false;
      }
      if (!q->tip) {
        for (j = 1; j <= n - 2; j++)
          putc('-', outfile);
        if (q->number - numsp >= 10)
          fprintf(outfile, "%2ld", q->number - numsp);
        else
          fprintf(outfile, "-%ld", q->number - numsp);
        extra = true;
      } else {
        for (j = 1; j < n; j++)
          putc('-', outfile);
      }
    } else if (!p->tip) {
      if (last->ycoord > i && first->ycoord < i && i != p->ycoord) {
        putc('!', outfile);
        for (j = 1; j < n; j++)
          putc(' ', outfile);
      } else {
        for (j = 1; j <= n; j++)
          putc(' ', outfile);
        trif = false;
      }
    }
    if (q != p)
      p = q;
  } while (!done);
  if (p->ycoord == i && p->tip) {
    /*
    for (j = 0; j < namelength; j++)
      putc(p->nayme[j], outfile);
      */
    fprintf(outfile, "%s", p->nayme);
  }
  putc('\n', outfile);
}  /* drawline */

void printree()
{
  /* prints out diagram of the tree */

  long tipy = 1;
  double scale, tipmax = 0.0;
  long i;

  if (!treeprint)
    return;
  putc('\n', outfile);
  coordinates(curtree.start, 0.0, &tipy,&tipmax);
  scale = 1.0 / (long)(tipmax + 1.000);
  for (i = 1; i <= tipy - down; i++)
    drawline(i, scale);
  putc('\n', outfile);
}  /* printree */


void describe(p)
node *p;
{
  /* print out information for one branch */
  long i;
  node *q;

  q = p->back;
  fprintf(outfile, "%4ld          ", q->number - numsp);
  if (p->tip) {
    /*
    for (i = 0; i < namelength; i++)
      putc(p->nayme[i], outfile);
      */
    fprintf(outfile, "%s", p->nayme);
  } else
    fprintf(outfile, "%4ld      ", p->number - numsp);
  fprintf(outfile, "%15.5f\n", q->v);
  if (!p->tip) {
    describe(p->next->back);
    describe(p->next->next->back);
  }
}  /* describe */

void summarize()
{
  /* print out branch lengths etc. */
  putc('\n', outfile);
  if (njoin) {
    fprintf(outfile, "remember:");
    if (outgropt)
      //fprintf(outfile, " (although rooted by outgroup)");
    fprintf(outfile, " this is an unrooted tree!\n");
  }
  fprintf(outfile, "\nBetween        And            Length\n");
  fprintf(outfile, "-------        ---            ------\n");
  describe(curtree.start->next->back);
  describe(curtree.start->next->next->back);
  if (njoin)
    describe(curtree.start->back);
  fprintf(outfile, "\n\n");
}  /* summarize */

void treeout(p)
node *p;
{
  /* write out file with representation of final tree */
  long i, n, w;
  Char c;
  double x;

  if (p->tip) {
    n = 0;
    /*
    for (i = 1; i <= namelength; i++) {
      if (p->nayme[i - 1] != ' ')
        n = i;
    }
    for (i = 0; i < n; i++) {
      c = p->nayme[i];
      if ( (c == ' ') || (c == '(') || (c == ')') )
        c = '_';
      putc(c, treefile);
    }
    */
    for (i = 1; i <= namelength; ++i) {
      if (p->nayme[i-1] != '\0')
	n = i;
    }
    fprintf(treefile,"%s", p->nayme);
    col += n;
  } else {
    putc('(', treefile);
    col++;
    treeout(p->next->back);
    putc(',', treefile);
    col++;
    if (col > 55) {
      putc('\n', treefile);
      col = 0;
    }
    treeout(p->next->next->back);
    if (p == curtree.start && njoin) {
      putc(',', treefile);
      treeout(p->back);
    }
    putc(')', treefile);
    col++;
  }
  x = p->v;
  if (x > 0.0)
    w = (long)(0.43429448222 * log(x));
  else if (x == 0.0)
    w = 0;
  else
    w = (long)(0.43429448222 * log(-x)) + 1;
  if (w < 0)
    w = 0;
  if (p == curtree.start)
    fprintf(treefile, ";\n");
  else {
    fprintf(treefile, ":%*.5f", (int)(w + 7), x);
    col += w + 8;
  }
}  /* treeout */

void nodelabel(isnode)
boolean isnode;
{
  if (isnode)
    printf("NODE");
  else
    printf("OTU ");
}  /* nodelabel */

Local Void jointree()
{
  /* calculate the tree */
  long nc, nextnode, mini, minj, i, j, ia, ja, ii, jj, nude, iter;
  double diq, djq, dij, fotu2, total, tmin, dio, djo, bi, bj, bk, dmin, da;
  long el[3];
  vector av;
  intvector oc;

  for (i = 0; i <= numsp - 2; i++) {
    for (j = i + 1; j < numsp; j++) {
      da = (x[i][j] + x[j][i]) / 2.0;
      x[i][j] = da;
      x[j][i] = da;
    }
  }
  /* First initialization */
  fotu2 = numsp - 2.0;
  nextnode = numsp + 1;
  av = (vector)Malloc(numsp*sizeof(double));
  oc = (intvector)Malloc(numsp*sizeof(long));
  for (i = 0; i < numsp; i++) {
    av[i] = 0.0;
    oc[i] = 1;
  }
  /* Enter the main cycle */
  if (njoin)
    iter = numsp - 3;
  else
    iter = numsp - 1;
  for (nc = 1; nc <= iter; nc++) {
    for (j = 2; j <= numsp; j++) {
      for (i = 0; i <= j - 2; i++)
        x[j - 1][i] = x[i][j - 1];
    }
    tmin = 99999.0;
    /* Compute sij and minimize */
    for (ja = 2; ja <= numsp; ja++) {
      jj = enterorder[ja - 1];
      if (cluster[jj - 1] != NULL) {
        for (ia = 0; ia <= ja - 2; ia++) {
          ii = enterorder[ia];
          if (cluster[ii - 1] != NULL) {
            if (njoin) {
              diq = 0.0;
              djq = 0.0;
              dij = x[ii - 1][jj - 1];
              for (i = 0; i < numsp; i++) {
                diq += x[i][ii - 1];
                djq += x[i][jj - 1];
              }
              total = fotu2 * dij - diq - djq;
            } else
              total = x[ii - 1][jj - 1];
            if (total < tmin) {
              tmin = total;
              mini = ii;
              minj = jj;
            }
          }
        }
      }
    }
    /* compute lengths and print */
    if (njoin) {
      dio = 0.0;
      djo = 0.0;
      for (i = 0; i < numsp; i++) {
        dio += x[i][mini - 1];
        djo += x[i][minj - 1];
      }
      dmin = x[mini - 1][minj - 1];
      dio = (dio - dmin) / fotu2;
      djo = (djo - dmin) / fotu2;
      bi = (dmin + dio - djo) * 0.5;
      bj = dmin - bi;
      bi -= av[mini - 1];
      bj -= av[minj - 1];
    } else {
      bi = x[mini - 1][minj - 1] / 2.0 - av[mini - 1];
      bj = x[mini - 1][minj - 1] / 2.0 - av[minj - 1];
      av[mini - 1] += bi;
    }
    if (progress) {
      printf("CYCLE %3ld: ", iter - nc + 1);
      if (njoin)
        nodelabel(av[mini - 1] > 0.0);
      else
        nodelabel(oc[mini - 1] > 1.0);
      printf("%3ld (%10.5f) JOINS ", mini, bi);
      if (njoin)
        nodelabel(av[minj - 1] > 0.0);
      else
        nodelabel(oc[minj - 1] > 1.0);
      printf("%3ld (%10.5f)\n", minj, bj);
    }
    hookup(curtree.nodep[nextnode - 1]->next, cluster[mini - 1]);
    hookup(curtree.nodep[nextnode - 1]->next->next, cluster[minj - 1]);
    cluster[mini - 1]->v = bi;
    cluster[minj - 1]->v = bj;
    cluster[mini - 1]->back->v = bi;
    cluster[minj - 1]->back->v = bj;
    cluster[mini - 1] = curtree.nodep[nextnode - 1];
    cluster[minj - 1] = NULL;
    nextnode++;
    if (njoin)
      av[mini - 1] = dmin * 0.5;
    /* re-initialization */
    fotu2 -= 1.0;
    for (j = 0; j < numsp; j++) {
      if (cluster[j] != NULL) {
        if (njoin) {
          da = (x[mini - 1][j] + x[minj - 1][j]) * 0.5;
          if (mini - j - 1 < 0)
            x[mini - 1][j] = da;
          if (mini - j - 1 > 0)
            x[j][mini - 1] = da;
        } else {
          da = x[mini - 1][j] * oc[mini - 1] + x[minj - 1][j] * oc[minj - 1];
          da /= oc[mini - 1] + oc[minj - 1];
          x[mini - 1][j] = da;
          x[j][mini - 1] = da;
        }
      }
    }
    for (j = 0; j < numsp; j++) {
      x[minj - 1][j] = 0.0;
      x[j][minj - 1] = 0.0;
    }
    oc[mini - 1] += oc[minj - 1];
  }
  /* the last cycle */
  nude = 1;
  for (i = 1; i <= numsp; i++) {
    if (cluster[i - 1] != NULL) {
      el[nude - 1] = i;
      nude++;
    }
  }
  if (!njoin) {
    curtree.start = cluster[el[0] - 1];
    return;
  }
  bi = (x[el[0] - 1][el[1] - 1] + x[el[0] - 1][el[2] - 1] - x[el[1] - 1]
        [el[2] - 1]) * 0.5;
  bj = x[el[0] - 1][el[1] - 1] - bi;
  bk = x[el[0] - 1][el[2] - 1] - bi;
  bi -= av[el[0] - 1];
  bj -= av[el[1] - 1];
  bk -= av[el[2] - 1];
  if (progress) {
    printf("LAST CYCLE:\n");
    putchar(' ');
    nodelabel(av[el[0] - 1] > 0.0);
    printf("%3ld  (%10.5f) JOINS ", el[0], bi);
    nodelabel(av[el[1] - 1] > 0.0);
    printf("%3ld  (%10.5f) JOINS ", el[1], bj);
    nodelabel(av[el[2] - 1] > 0.0);
    printf("%3ld  (%10.5f)\n", el[2], bk);
  }
  hookup(curtree.nodep[nextnode - 1], cluster[el[0] - 1]);
  hookup(curtree.nodep[nextnode - 1]->next, cluster[el[1] - 1]);
  hookup(curtree.nodep[nextnode - 1]->next->next, cluster[el[2] - 1]);
  cluster[el[0] - 1]->v = bi;
  cluster[el[1] - 1]->v = bj;
  cluster[el[2] - 1]->v = bk;
  cluster[el[0] - 1]->back->v = bi;
  cluster[el[1] - 1]->back->v = bj;
  cluster[el[2] - 1]->back->v = bk;
  curtree.start = cluster[el[0] - 1]->back;
  free(av);
  free(oc);
}  /* jointree */


void maketree()
{
  /* construct the tree */
  long i, j, k;

  getdata();
  if (progress)
    putchar('\n');
  if (ith == 1)
    setuptree(&curtree);
  for (i = 1; i <= numsp; i++)
    enterorder[i - 1] = i;
  if (jumble) {
    for (i = 0; i < numsp; i++) {
      j = (long)(randum(seed) * numsp) + 1;
      k = enterorder[j - 1];
      enterorder[j - 1] = enterorder[i];
      enterorder[i] = k;
    }
  }
  for (i = 0; i < numsp; i++)
    cluster[i] = curtree.nodep[i];
  jointree();
  if (outgropt)
    curtree.start = curtree.nodep[outgrno - 1]->back;
  printree();
  if (treeprint)
    summarize();
  if (trout) {
    col = 0;
    treeout(curtree.start);
  }
  if (progress) {
    printf("\nUPGMA Tree written on aniupgma\.txt file\n\n");
    if (trout)
      printf("Cluster written on treefile file\n\n");
  }
}  /* maketree */


main(argc, argv)
int argc;
Char *argv[];
{  /* main program */
char infilename[100],outfilename[100],trfilename[100];
#ifdef MAC
  macsetup("Neighbor","");
  argv[0] = "Neighbor";
#endif
  openfile(&infile,"anir.txt","r",argv[0],infilename);
  openfile(&outfile,"aniupgma.txt","w",argv[0],outfilename);
  ibmpc = ibmpc0;
  ansi = ansi0;
  vt52 = vt520;
  mulsets = false;
  datasets = 1;
  doinit();
  if (trout)
    openfile(&treefile,TREEFILE,"w",argv[0],trfilename);
  x = (vector *)Malloc(numsp*sizeof(vector));
  for (j = 0; j < numsp; j++)
    x[j] = (vector)Malloc(numsp*sizeof(double));
  reps = (intvector *)Malloc(numsp*sizeof(intvector));
  for (j = 0; j < numsp; j++)
    reps[j] = (intvector)Malloc(numsp*sizeof(long));
  names = (naym *)Malloc(numsp*sizeof(naym));
  enterorder = (long *)Malloc(numsp*sizeof(long));
  cluster = (node **)Malloc(numsp*sizeof(node *));
  ith = 1;
  while (ith <= datasets) {
    if (datasets > 1) {
      fprintf(outfile, "Data set # %ld:\n",ith);
      if (progress)
        printf("Data set # %ld:\n",ith);
    }
    getinput();
    maketree();
    if (eoln(infile)) {
      fscanf(infile, "%*[^\n]");
      getc(infile);
    }
    ith++;
  }
  FClose(infile);
  FClose(outfile);
  FClose(treefile);
#ifdef MAC
  fixmacfile(outfilename);
  fixmacfile(trfilename);
#endif
  exit(0);
}

int eof(f)
FILE *f;
{
    register int ch;

    if (feof(f))
        return 1;
    if (f == stdin)
        return 0;
    ch = getc(f);
    if (ch == EOF)
        return 1;
    ungetc(ch, f);
    return 0;
}


int eoln(f)
FILE *f;
{
    register int ch;

    ch = getc(f);
    if (ch == EOF)
        return 1;
    ungetc(ch, f);
    return (ch == '\n');
}

void memerror()
{
printf("Error allocating memory\n");
exit(-1);
}

MALLOCRETURN *mymalloc(x)
long x;
{
MALLOCRETURN *mem;
mem = (MALLOCRETURN *)malloc(x);
if (!mem)
  memerror();
else
  return (MALLOCRETURN *)mem;
}

